function [output,SortResults,OutTable] = estimate_k(k_target,pars,Nu,M_full,W_full,random_state,verbose,xi_grid)     
    
target_ratio                            = 1 + pars.theta/(pars.r + k_target);
if nargin == 8
    Ngrid                                   = numel(xi_grid);
else
    Ngrid                                   = 1e1+1;
    xi_grid                                 = linspace(0.00700,0.01300,Ngrid)';
end
sigma_grid                              = xi_grid*target_ratio;
SortResults                             = NaN(Ngrid,8);
all_results                             = cell(Ngrid,1);
Nu_initial                              = Nu;
pars_initial                            = pars;

for n=1:Ngrid
    if verbose == 1
        fprintf('Run # %.0f \n',n)
        fprinct('********** \n');
    end
    setup_inputs.k                          = k_target;
    setup_inputs.xi                         = xi_grid(n);
    setup_inputs.sigma                      = sigma_grid(n);
    setup_inputs.C                          = (1/2)*( Nu_initial.lb_global(5) + min(2*sqrt(Nu_initial.N/pars_initial.theta)*setup_inputs.xi,Nu_initial.ub_global(5)) );
    setup_inputs.M_e                        = (1/2)*( max(Nu_initial.lb_global(4),pars_initial.M_c - sqrt(Nu_initial.N/pars_initial.theta)*setup_inputs.xi) + min(Nu_initial.ub_global(4),pars_initial.M_c + sqrt(Nu_initial.N/pars_initial.theta)*setup_inputs.xi-setup_inputs.C) );
    setup_inputs.S                          = 0.246;
    if strcmp(Nu.which_moments,'baseline')
    [pars,Nu]                               = setup_function(setup_inputs,random_state,'baseline');
    else
    [pars,Nu]                               = setup_function(setup_inputs,random_state,'');
    end
    all_results{n}                          = estimate(pars,M_full,W_full,Nu,0);
    SortResults(n,1)                        = n;
    SortResults(n,2)                        = all_results{n}.ObjFinal;
    SortResults(n,3:end)                    = all_results{n}.XiFinalAllPars;
    all_results{n}.setup_inputs             = pars;
    all_results{n}.pars                     = pars;
end

SortResults                             = sortrows(SortResults,2);
output                                  = all_results;
OutTable                                = table( ...
                                                strvcat(num2str(SortResults(1,5),'%.5f'),[' ']),...
                                                strvcat(num2str(SortResults(1,2),'%.3f'),[' ']),...
                                                strvcat(num2str(SortResults(1,3),'%.5f'),[' ']),...
                                                strvcat(num2str(SortResults(1,4),'%.5f'),[' ']),...
                                                strvcat(num2str(SortResults(1,6),'%.5f'),[' ']),...
                                                strvcat(num2str(SortResults(1,7),'%.5f'),[' ']),...
                                                strvcat(num2str(SortResults(1,8),'%.5f'),[' ']),...
                                                [SortResults(1,1);0]);
OutTable                                = OutTable(1,:);
OutTable.Properties.VariableNames       = ["k" "Obj." "xi" "sigma" "S" "M_e" "C" "Run #" ];

if verbose
    disp(OutTable)
end

end